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Abstract 

We present the calculations of the complete NLO inclusive total cross sections for pair production 
of neutral Higgs bosons through hb annihilation in the minimal supersymmetric standard model 
at the CERN Large Hadron Collider. In our calculations, we used both the DREG scheme and 
the DRED scheme and found that the NLO total cross sections in these two schemes are the 
same. Our results show that the 66-annihilation contributions can exceed those of gg fusion and qq 
annihilation for h^H^, A^h^ and A^H^ productions when tan/3 is large. In the case of /i > 0, the 
NLO corrections enhance the LO total cross sections significantly, reaching a few tens of percent, 
while for /i < 0, the corrections are relatively small, and are negative in most of parameter space. 
Moreover, the NLO QCD corrections reduce the dependence of the total cross sections on the 
renormalization/factorization scale, especially for ^ < 0. We also used the CTEQ6.1 PDF sets 
to estimate the uncertainty of LO and NLO total cross sections, and found that the uncertainty 
arising from the choice of PDFs increases with the increasing m^o . 

PACS numbers: 12.60.Jv, 12.38.Bx, 13.85.Fb 
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I. INTRODUCTION 



The Higgs mechanism plays a key role for spontaneous breaking of the electroweak symme- 
try both in the standard model (SM) and in the minimal supersymmetric (MSSM) extension 
of the SM ^. Therefore, the search for Higgs bosons becomes one of the prime tasks in 
future high-energy experiments, especially at the CERN Large Hadron Collider (LHC), with 
a/S* = 14 TeV and a luminosity of 100 fb^^ per year In the SM, only one Higgs doublet 
is introduced, and the neutral CP-even Higgs boson mass is basically a free parameter with 
a theoretical upper bound of rriH < 600 - 800 GeV [3] and a LEP2 experimental lower 
bound of rriH > 114.4 GeV |4:j. In the MSSM, two Higgs doublets are required in order to 
preserve supersymmetry (SUSY), and consequently the model predicts five physical Higgs 
bosons: the neutral CP-even ones and H^, the neutral CP-odd one A^, and the charged 
ones H^. The h^, which behaves like the SM one in the decoupling region (m^o ^ rnzo), 
is the lightest, and its mass is constrained by a theoretical upper bound of rriho < 140 GeV 
when including the radiative corrections ^. The analyses in 1^ indicate that the /i^ boson 
can not escape detection at the LHC, and that in large areas of the parameter space, more 
than one Higgs particle in the MSSM can possibly be found, which is an exciting result, 
since the discovery of any additional Higgs bosons will be direct evidence of physics beyond 
the SM. 

At the LHC, a neutral Higgs boson can be produced through followingmechanisms: 
gluon fusion ^ weak boson fusion qq — >• qgV*V* — >• qqh^/qqH^ associated 
production with weak bosons ^j, pair production [3, llll IlE 13 > ^ind associated production 
with a tt pair gg/qq ^ tt(j) In the MSSM, since the couplings between Higgs bosons and 
b quarks can be enhanced by large values of tan (3, the ratio of the vacuum expectation values 
of the two Higgs doublets, Higgs bosons will also be copiously produced in association with b 
quarks at the LHC. Except for qq —>■ bbh^, the other relevant production mechanisms depend 
on the final state being observed il5]. For inclusive Higgs production, the lowest order 
process is bb [3, 

and the convergence of the perturbative expansion is improved by 
summing the coUinear logarithms to all orders through the use of b quark parton distributions 
with an appropriate factorization scale. However, if at least one high-p^ b quark is required 
to be observed, the leading partonic process is gb — > bh^ and if two high-p^ b quarks 
are required, the leading subprocess is gg — * bbh^ 
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Studying the pair production of neutral Higgs bosons may be an important way to probe 
the trihnear neutral Higgs boson couplings, which can distinguish between the SM and the 
MSSM. In the SM, Higgs boson pair production is dominated by gg fusion mediated via 
heavy-quark loops, while the contribution of qq annihilation is greatly suppressed by the 
absence of the HHZ coupling and the smallness of the Hqq (g = u,d,s, c, h) couplings. 
In the MSSM, gg fusion for the pair production of neutral Higgs bosons can be mediated 
via both quark loops 11 1 and squark loops and the existence of h^A^Z and 



H^A^Z couplings at the tree level leads to h^A^ and H^A^ associated productions through qq 
annihilations (Drell-Yan-like processes) . Moreover, since the (pbb couplings can be greatly 
enhanced by large values of tan /?, there are potentially important contributions arising from 
hb annihilation to pair production of neutral Higgs bosons, which have been studied at the 
leading-order (LO) [13]. However, the LO predictions generally have a large uncertainty 
due to scale and PDF choices. In this paper, we present the complete next-to-leading order 
(NLO) QCD (including SUSY-QCD) calculation for the cross sections for pair production 
of neutral Higgs boson through hh annihilation at the LHC Similar to single Higgs boson 
production, for the inclusive production the use of h quark parton distributions at the LO 
will improve the convergence of the perturbative expansion. For simplicity, we neglect the 
bottom quark mass except in the Yukawa couplings, which is valid in all diagrams where the 
bottom quark is an initial state parton, according to the simplified Aivazis-CoUins-Olness- 
Tung (ACOT) scheme For regularization of the ultraviolet (UV), soft and coUinear 

divergences, both the dimensional regularization (DREG) approach j2o| (with naive 75 [2^) 
and the dimensional reduction (DRED) scheme 3| are used in our calculations providing 
a cross check. 

This paper is organized as follows. In Sect. II we show the analytic results for the LO 
cross sections proceeding through hb annihilation. In Sect. HI we present the details of the 
calculations of both the virtual and real parts of the NLO QCD corrections, and compare the 
results using DREG with those using DRED. In Sect. IV we give the numerical predictions 
for inclusive and differential cross sections at the LHC. The relevant coupling constants and 
the lengthy analytic expressions are summarized in Appendices A, B and C. 
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II. LEADING ORDER PAIR PRODUCTION OF NEUTRAL HIGGS BOSONS 



The tree-level Feynman diagrams for the subprocess b{pa)b{pb) — > Hi{pi)Hj{p2), where 
Hi=i,2,3 = 1 hP-iA^, are shown in Fig. ^ and its LO amplitude in n = 4 — 2e dimensions is 

= ^^\Ml^ + + Mlfl (2.1) 



with 



M-j' = -2^ — v{pb){akPL + a^PuMPa) 



t 



M>/ = — - — v{pb){a.ja*PL + aia*PR) ^iu{pa), 



M^f = 'l^v{pb){a,a]PL + a.ja*PR) f^MPa), (2.2) 

where = G°, cw = cos9w, Yb = m^/ {\/2mw cos (3), Pl,r = (1 T 75)/2, C^l = 
—1/2 + sin^^v^/S, CbR = siTL^Ow/S, and fir is a mass parameter introduced to keep the 
coupling constant g dimensionless. sx = s — m\ + ifnx^x is the denominator of the prop- 
agator of particle X with mass nix and total decay width Fx- Ckij, au and denote the 
coefficients appearing in the HkHiHj, Hkbb and Z^A^Hj couplings, respectively, and their 
explicit expressions are shown in Appendix A. Mandelstam variables s, t and u are defined 
as follows 

S = {Pa+ Pbf, t={pa- Plf, U= {pa- P2f ■ (2.3) 

The above amplitude and all of the other calculations in this paper are carried out in 
t'Hooft-Feynman gauge. 

After the n-dimensional phase space integration, the LO parton level differential cross 
sections are 

d'^^fj 1 7r5e ,t'u' - sm]j , 2 ^r^\ I ^ 



x5is + t + u-m%-m%)\M,^\' (2.4) 

where 5"^ = (47r)~^+% t' = t — m'jj_, u' = u — mj^ , the factor 1/{1 + Sij) accounts for identical- 
particle symmetrization when Hi = Hj. \Mlj\ is the LO amplitude squared, where the 
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colors and spins of the out going particles have been summed over, and the colors and spins 



of the incoming ones have been averaged over. The explicit expression for \M^\ is 



^(-W k,l=l ^Hk^H, ^ 



+ 5, 



i3 



{Zj ) (C^^ + C^^) 2iY'b CbRaja* + CbLO'iO'*j CbL^ja* + CbRaia* , 



t 



u 



[2.5) 



The LO total cross section at the LHC is obtained by convoluting the parton level cross 
section with the parton distribution functions (PDFs) Gb^b/p the proton: 



a, 



J dxidx2[Gb/pixi, Hf)G-b/p{x2, fJ^f) + (xi ^ Xs)]^^ 



(2.6) 



where /x/ is the factorization scale. 



III. NEXT-TO-LEADING ORDER CALCULATIONS 

The NLO corrections to pair production of neutral Higgs bosons through bb annihilation 
consist of the virtual corrections, generated by loop diagrams of colored particles, and the 
real corrections with the radiation of a real gluon or a massless (anti)bottom quark. For 
both virtual and real corrections, we will first present the results in the DREG scheme, and 
then in the DRED scheme and compare them. 



A. virtual corrections 

The Feynman diagrams for the virtual corrections to bb HiHj are shown in Fig. |21 In 
order to remove the UV divergences, we renormalize the bottom quark mass in the Yukawa 
couplin gs a nd the wave function of the bottom quark, adopting the on-shell renormalization 
scheme j22|. The relations between the bare bottom quark mass mbo, the bare wave function 
ipbo and their relevant renormalization constants 6mb, 5ZbL(R) are defined as 

rribo = mb + Srub, 

iJbo = (1 + 6ZbL)hbL + (1 + 6ZbR)hbR. (3.1) 
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Calculating the self-energy diagrams in Fig. |21 we obtain the explicit expressions for Smj, 
and 5Z^L{Ry 



^ = -^cJsB.iml 0,ml)-2 + Y: \b, - ^ sin 2^,(-l)^5o 

1=1 
i=i 

where Cp = 4/3, i?o,i are the two-point integrals 13], ^i^^ are the sbottom masses, is 
the gluino mass, and i?^ is a 2 x 2 matrix defined to rotate the sbottom current eigenstates 
into the mass eigenstates: 

(3.2) 



V 



R' 



\ 



>R 



COS 6i sin 01 



\ 



sin cos 9i 



with < ^^f, < 77 by convention. Correspondingly, the mass eigenvalues m^^^ and (with 



6i < "^bj are given by 



"^'^^ ° I = R''Ml{Ry, Ml 
m?^ 



/ 2 



(3.3) 



with 



mi 
ml 

bR 



M% + ml + ml cos 2(3Ci 



bL, 



Ml + ml -ml cos2pCbR, 



ab = Ab — /itan/3. 



(3.4) 



Here M^ is the sbottom mass matrix. Mq jj and Ab are soft SUSY-breaking parameters and 
/i is the higgsino mass parameter. 

The renormalized virtual amplitudes can be written as 



(3.5) 



Here M^^"^^" contains the self-energy, vertex and box corrections, and can be written as 



j^unren ^ ^ ^^2^2 -(^^^ [^.p^ ^ ^ap^^ ^^(^ap^ ^ /rP«)]«(p.) 



(3.6) 
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where a denotes the corresponding diagram in Fig. |21 and (/ = 1,2,3,4) are the form 
factors given exphcitly in Appendix B. M™" is the corresponding counterterm, and can be 

respectively: 



separated into M^°^'^^\ M^""'^^^ and M^°"'^^\ i.e. the counterterms for s, t and u channels, 



with 



M^P = M^p^'^ + M^°"(*) + Mf;"(") (3.7) 



Zcw rrib / f^^^ Sh^, 



^coniu) ^ ^^[a,a*(2^ + 5ZbL)PL + «^«K2^ + 5ZbR)PR]u{pa). (3.8) 

The 0{as) virtual corrections to the differential cross section can be expressed as 



x5{s + t + u-m%-m%)2 Re(Mp^, (3.9) 

where the renormalized amplitude is UV finite, but still contains the infrared (IR) 
divergences, and is given by 



with 

AY = -Cp, AY = -^Cf. (3.11) 

The coefficients and A^ are constants, and similar to those in the pure Drell-Yan-like 
processes (without color particles in the final states). These IR divergences include the 
soft divergences and the collinear divergences. The soft divergences will be cancelled after 
adding the real corrections, and the remaining collinear divergences can be absorbed into 
the redefinition of PDFs j25|, which will be discussed in the following subsections. 

When recalculating the above virtual corrections in the DRED scheme, one finds that 
6ZbL and SZbR remain unchanged, however, 6mb and the form factors have shifts which are. 



respectively, given by 



Smh\ f 5m, 



mi, J DREG \ mb J DRED 4:71 



a 



+ -rCF, (3.12) 



and 



Jl DREG ~ Jl DRED ' 

a=a a=a k=l '^W^Hf^ 

J 2 DREG — 2^ J2 DRED ' 

a=a a=a k=l ^W^H^, 



fs DREG ~ fs DRED 5 / ^2 \ i.'^'^i'^j 

a=a a=a '^WK'^ '"'Z) ^" 

d-^iZ^ ChR 3Y'^ 

X! ^4 DijBG — DRED ~ „2 /' „ ™2 N ——{uaja^ — taia^). (3.13) 



Thus it is easy to obtain the following relations: 

^^DREG = M^DRED " ^^^^l (3.14) 

""Idreg = ""Idred - ^^^""S + ^("')' (3-15) 
where and af^ are independent of the choice of schemes. 

B. Real gluon emission 

The feynman diagrams for the real gluon emission process b{pa)b{pb) — > Hi{pi)Hj{p2) + 
g{p^) are shown in Fig. |H| 

The phase space integration for the real gluon emission will produce soft and collinear 
singularities, which can be conveniently isolated by slicing the phase space into different 
regions using suitable cut-offs. In this paper, we use the two cut-off phase space slicing 



method [26|, which introduces two arbitrary small cut-offs, i.e. soft cut-off 5s and collinear 
one (5c, to decompose the three-body phase space into three regions. 

First, the phase space is separated into two regions by the soft cut-off 5^, according to 
whether the energy of the emitted gluon is soft, i.e. < 6sy/s/2, or hard, i.e. E^, > 5s\fs/2. 
Correspondingly, the parton level real cross section can be written as 

where af^ and afj are the contributions from the soft and hard regions, respectively, 
contains all the soft divergences, which can explicitly be obtained after the integration over 
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the phase space of the emitted gluon. Next, in order to isolate the remaining colhnear 
divergences from a^, the colhnear cut-off 5c is introduced to further split the hard gluon 
phase space into two regions, according to whether the Mandelstam variables Mi 2 = {Pa,b — 
p^Y satisfy the colhnear condition —5cS < Ui^2 < or not. We then have 

^ij = ^ij ' (3.17) 

where the hard colhnear part af^'-^ contains the colhnear divergences, which also can explicitly 
be obtained after the integration over the phase space of the emitted gluon. And the hard 
non-coUinear part o"^*" is finite and can be numerically computed using standard Monte- 
Carlo integration techniques and can be written in the form 

2s' 



daJ^ = -|M^| dV,. (3.18) 



Here dT^ is the hard non-collinear region of the three-body phase space, and the explicit 



expressions for {M-'jl are given in Appendix C. 

In the next two subsections, we will discuss in detail the soft and hard colhnear gluon 
emission. 

1. Soft gluon emission 

In the limit that the energy of the emitted gluon becomes small, i.e. E3 < 6s\^/2, the 



amplitude squared \M^\ can simply be factorized into the Born amplitude squared times 
an eikonal factor ^eik'- 



iM^ibb ^ H,Hj + g)f {A7las^^l')\M,^f<^e^k, (3.19) 



where the eikonal factor $ejfc is given by 



<^eik = Cf- -. (3.20) 

(Pa ■P3)(P6 -Ps) 

Moreover, the phase space in the soft limit can also be factorized: 

dT^{hh HiHj + g) dT2{hh HiHj)dS, (3.21) 
Here dS is the integration over the phase space of the soft gluon, and is given by 

= 2(2^0^ i dE,Er^dn,_2.. (3.22) 
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The parton level cross section in the soft region can then be expressed as 

= {Anasfil') J dT^W^^ J dS<l>eik. (3.23) 
after the integration over the soft gluon phase space, 



Using the approach of Ref. 
Eq. ()3.23|) becomes 



Ois r(l — e) / 47r/i 



27rr('l - 2e) 



4* 4'' 



(3.24) 



with 



Al = 2Cf , Al = -4Cf In 5,, = ACp In^ 5, 



(3.25) 



These coefficients are the same as the ones in pure Drell-Yan-like processes, as expected. 



2. Hard coUinear gluon emission 

In the hard collinear region, E3 > 6s\^/2 and —6cS < ui^2 < 0, the emitted hard gluon 
is collinear to one of the incoming partons. As a consequence of the factorization theorems 
[2^, the amplitude squared for bb HiHj + g can be factorized into the product of the 
Born amplitude squared and the Altarelli-Parisi splitting function for b{b) b(b)g 29 



iM^ibb ^ H,H^ + g)\ 



■2 colline 



(47ra,/i2^)|M^| 



■2/-2Pfeb(z,e) , -2Pii{z,e 



+ 



3- 

(3.26) 



ZUi ZU2 

Here z denotes the fraction of incoming parton b{b)^s momentum carried by parton b{b) 
with the emitted gluon taking a fraction (1 — z), and Pij{z, e) are the unregulated splitting 
functions in n = 4 — 2e dimensions for < z < 1, which can be related to the usual 
Altarelli-Parisi sphtting kernels j^l as follows: Pij{z,e) = Pij{z) + ePlj{z), explicitly 



P,f,{z) = Pii{z 



Cp-^ + CfU{1-z), 
1 — z 2 



PLiz) = Pij,{z) = -C^(l -z)+ Cf-5{1 - z). 



(3.27) 
(3.28) 



Moreover, the three-body phase space can also be factorized in the collinear limit, and, for 
example, in the limit —5cS < Mi < it has the following form 

(4vr) 



dV^{bb H,Hj + g) dV2{bb H,Hj- s' = zs) 



levr^rci - e 



-dzdui[-{l - z)ui]X^.29) 
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Here the two-body phase space is evaluated at a squared parton-parton energy of zs. Thus 
the three-body cross section in the hard colhnear region is given by 

'as r(l — e) AiT^l 



da, 



HC 



where G 



6(5)/pl 



27rr(l-2e)' s 

+Pll{z,e)Gii.p{xi/z)Gi,ip{x2) + {xi ^ x^ 

x\ is the bare PDF. 



Pbb{z,e)Gb/p{xi/z)Gi/p{x2 
dz , \ — z 



) ''dxidx2j (3.30) 



C. Massless 6(6) emission 

In addition to real gluon emission, a second set of real emission corrections to the inclusive 
cross section for pp — * HiHj at NLO involves the processes with an additional massless h{h) 
in the final state: 



HiHjb, 



gb HiHjb. 



The relevant feynman diagrams for massless b emission are shown in Fig. 01 and the diagrams 
for 6-emission are similar and omitted here. 

Since the contributions from real massless b(b) emission contain initial state collinear 

n 

singularities, we also need to use the two cut-off phase space slicing method j26|] to isolate 
these collinear divergences. But we only split the phase space into two regions, because there 
are no soft divergences. Consequently, using the approach in Ref. the cross sections for 
the processes with an additional massless b(b) in the final state can be expressed as 



datf'' 



crfj{ga HiHj + X)[Gg/p{xi)Ga/p{x2) + {xi ^ X2)]dxidx2 

a=b,b 



+da. 



B 



as r(l - e) Anfi: 



2 ■ 

r \e 



A), 

e 



.27rr(l-2e)' s 
+Gb/p{xi)Pig{z,e)Gg/p{x2/z) + (xi 



Pbg{z, t)Gglp{Xi/ z)Gb/p{x2) 



X2 



dz / \ — z 
z \ z 



dxidx2, 



(3.31) 



where 



\z' + (1 - z?] 



Hg{z)=Pig{z) = -z{l-z). (3.32) 



The first term in Eq. ()3.31|) represents the non-coUinear cross sections for the two processes, 
which can also be written in the form {a = b, b) 

72 



da, 



1 
Ys 



(3.33) 
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where dT^ is the three body phase space in the non-colhriear region. The exphcit expressions 

2 ^2 

for |Mj^"| can be obtained from {M^jl by crossing symmetry. The second term in Eq. ()3.31|) 
represents the colhnear singular cross sections. 



D. Mass factorization 



As mentioned above, after adding the renormahzed virtual corrections and the real cor- 
rections, the parton level cross sections still contain colhnear divergences, which can be 
absorbed into the redefinition of the PDFs at NLO, in general called mass factorization j^]. 
This procedure in practice means that first we convolute the partonic cross section with the 
bare PDF Ga/p{x), and then use the renormahzed PDF Ga/p{x, fij) to replace Ga/p{x). In 
the MS convention, the scale dependent PDF Ga/p{x, ^f) is given by 



Ga/pix,lXf) = Ga/p{x) + J2{ ) 

/3 ^ 



«s r(l — e) / iiTfil 



1 dz 



Pc.p{z)Gp,p{x/z). (3.34) 



.27rr(l - 2e) V ^i) 

This replacement will produce a colhnear singular counterterm, which is combined with the 
hard colhnear contributions to result in, as the definition in Ref. 12^], the 0{as) expression 
for the remaining colhnear contribution: 

'as r(l - e) Mvr/i^^ 



da. 



coll 



daf^ 



+ E 

a=b,b 



27rr(l-2e)V s 

A'^^{a ag) 



{Gb/p{xi, fif)Gf^/.p{x2, fif) + Gh/p{xi, fif)Gf,/p{x2, fif) 



+ A'^'ia ag) 



Gb/p{xi, fif)Gf,/p{x2, fif) 



+ (xi ^ X2)}dXidX2, 



(3.35) 



where 



with 



ATib ^ bg) = ATib bg) = C^(21n5, + 3/2), 
A- = Arin(4), 

Ga/p{x,flf) = XI / " —Gl3/p{x/y,^f)Pa(5{y) 

R -^^ y 



PaM = Pap in(5e^4) - Kp{y)- 
y /^/ 



(3.36) 
(3.37) 

(3.38) 



(3.39) 
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Finally, the NLO total cross section for pp — * HiHj in the MS factorization scheme is 

= / {dxidx2 Gk/pixu fif)Gi/pix2, fif) + (xi ^ X2) (af. + al + af- + af^) + a'rf} 
+ J2 J dxidx2 Gg/p{xi, iJ,f)Ga/pix2, fif) + (xi ^ X2) o-^jiga HiHj + X). (3.40) 

a=b,b 

Note that the above expression contains no singularities since 2A2 + A| = and 2AY + AI + 
Af{b bg) + (6 bg) = 0. 

E. Real corrections and NLO total cross sections in the DRED scheme 

Above we gave the real corrections and NLO total cross sections in the DREG scheme, 
and next we show the corresponding results in the DRED scheme, where the contributions 
from soft gluon emission remain the same, while those from hard coUinear gluon emission 
and massless (anti) quark emission are different. These differences arise from the splitting 
functions and the PDFs. 

First, the splitting functions in the DRED scheme have no e parts, and we have 

Pij{z,e)DRED = Pij{z). (3.41) 

Then from Eq. ()3.35|) and ()3.41|) we obtain 

(y r^~^s^bf3 dv 

<'dreg = <DRED - / -Gp,p{x^/y, ^f)P\^{y)Gi,p{x2, /x/) 



+ 22 — Gp/p{x2/y,fJ^f)P'-bf3{y)Gb/p{xi,fif) + {xi<-^X2)}adxidx2, (3.42) 

p •'^2 y 

Secondly, the PDFs in the DRED and DREG schemes are related j^: 

Ga/p{x, fIf)DREG = Ga/p{x, ^f)DRED + P'af3{x/y). (3.43) 

Ztt ^ Jx y 

Substituting into the formula for the Born cross section, we obtain an additional difference 
at the 0{as) level arising from the PDFs: 

^kiDREG ^ ^iiDRED + of "lXI / -ZT^ P /P^^^l y ^ ^^ f) ^REdP' bp{y)Gljp{x2, ^^f) DRED 

ZTT p J XI y 

^ f-\ (^y 

+ J2 Gf3/p{x2/y, fif)DREDP'bp{y)Gb/p{xi, IJ,f)DRED + (xi ^ X2)}a^ dxidx2. 

p ''^2 y 

(3.44) 
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Equations ()3.42|) and ()3.44p are very similar except for the limits of the integral over y in the 
two expressions. Substituting the equations ()3.42|) . ()3.44|) and ()3.15|) into ()3.40p . we obtain 
the relation of the NLO total cross sections in the two schemes: 



*J DREG DRED + ^^Y. J^_^ ^ l3/p{Xl/y , l^f) P' bpiy)Gl/p{x2, ^f) 



p Ji-SsSbp y 

dy. 

p J^~^shi3 y 



+ ^ / — Gf3/p{x2/y, iJ,f)P'bpiy)Gi,/pixi, Hf) + (xi ^ X2)}crfjdxidx2 



a. 



-^C^aS + 0{al). (3.45) 
Using the explicit expressions for the e parts of the splitting functions P' , we find 

*J DREG V DRED^^ sJ- \^<J.-tyjj 

Therefore, the NLO total cross sections in the two schemes are the same. 



F. Differential cross section in transverse momentum and invariant mass 

In this subsection we present the differential cross section in the transverse momentum pt 
and the invariant mass. Using the notations defined in Ref . , the differential distribution 
with respect to pt and y of Hi for the processes 

p{Pa)+p{Pb) ^ H,{pi) + H,{p2)[+gip3)/Kp3)/Kp3)] (3.47) 

is given by 

- — J- = 2ptS J2 J _dxi J dx2XiGa/p{xi, fif)x2G/3/p{x2, /"/)^^7^' (3-48) 



where is the total center-of-mass energy of the collider, and 

2 T2U2 2 1 1 N 

Pt = ^ ^Hi^ ?/=oln(— ), 



2 



XI = X2 = ^ (3.49) 

S + U2 X1S+T2 

with T2 = {ph — PiY — Tn?Hi ^^'i ^2 = {Pa — PiY ~^ ''^Hi- '^^^ limits of integration over y and 
Pt are 

- y^'^'iPT) <y< 2/"^'^^(pr), 0<pT< (3.50) 
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with 



y'^"'''{pT) = arccoshf 



m 



Hi 



V2 + 



m 



= ^^{S + - m%y - Am%S . (3.51) 
The differential distribution with respect to the invariant mass MH^Hj is given by 

da 



S „ dr 



where dCf^.j^./dr is the parton luminosity: 



"'*-HiHj 



1 dx 

X 



with 



dr 



ill 



(3.52) 



(3.53) 



(3.54) 



IV. NUMERICAL RESULTS AND CONCLUSIONS 

In this section, we present the numerical results for total and differential cross sections 
for pair production of neutral Higgs bosons at the LHC In our numerical calculations, the 
SM parameters were taken to be aew('^H/) = 1/128, m^i/ = 80.425 GeV, mz = 91.1876 GeV 



and mt = 178.1 GeV 
and CTEQ6M PDFs 



. We used the two-loop evaluation for as{Q) '33\ {as{Mz) = 0.118) 
34| throughout the calculations of the NLO (LO) cross sections unless 
specified . Moreover, in order to improve the perturbative calculations, we took the running 
mass mb{Q) evaluated by the NLO formula j35[: 



mb{Q) = Uti{Q,mt)U5{mt,mf,)mb{'mb) 
with mb^nib) = 4.25 GeV 0|. The evolution factor Uf is 



(4.1) 



tts(<5i 



^ ^ asiQi) -as{Q2) j(f) 
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33-2/' 



jU) 



Air 

8982 - 504/ + 40/2 
3(33 -2/)2 



(4.2) 
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In addition, to also improve the perturbation calculations, especially for large tan/3, we 
made the following SUSY replacements in the tree-level couplings 35^ : 

rribiQ) 



TTlbiQ) 



Arrib 



2a, 



1 + Arrib ' 



(4.3) 



oo 2 2 2^ Mt tan plim^^ , m^^ , /i) 
SZTT'^myy sm p 



167r2 



i=l 



/xM2tan/3^ (4)2/(m,-^, M2, /i) + M2, /i) 



(4.4) 



with 



I(a,b, c) 



1 



—^{a'b' log ^ + fo^c^ log ^ + c^a^ log 



(4.5) 



(a2-fe2)(52_c2)(a2 

It is necessary, to avoid double counting, to subtract these (SUSY-)QCD corrections from 
the renormalization constant 6mh in the following numerical calculations. 

For the MSSM parameters, we chose mi, mo, Aq, tan/5 and the sign of /i as input 
parameters, where mi, mo and Aq are the universal gaugino mass , scalar mass at the GUT 
scale and the trilinear soft breaking parameter in the superpotential terms, respectively. 
Specifically, we took mi = 170 GeV, Aq = 200 GeV, and tuned mo to obtain the desired 
value of rriy^o, while tan/5 and sign of /i are free. All other MSSM parameters are determined 
in the minimal supergravity (mSUGRA) scenario by the program package SUSPECT 2.3 
js^l- In particular, we used running DR Higgs masses at m^, which include the full one-loop 
corrections, as well as the two-loop corrections controlle d by the strong gauge coupling and 
the Yukawa couplings of the third generation fermions 38[. Moreover, since an s-channel 
resonance can occur in the process — >■ h^h^ + X when m/^o > 2m/jO, we adopted the 



program package HDECAY 3.101 39[ to determine the total decay widths of the Higgs 
bosons. For example, in the case of tan/5 = 40, m^o = 250 GeV and /i < 0, we have 
m^o = 108.7 GeV, m^o = 250.3 GeV, F/^o = 8 MeV, Fj^o = 13.1 GeV, F^o = 13.5 GeV, 
m-g = 454 GeV, m^^ = 422 GeV, and m^^ = 493 GeV. 

For the renormalization and factorization scales, we always chose Hr = m^v = {rriHi + 
mH^)/2 and /i/ = mav unless specified. 

In Fig. |3[ we chose A^H^ production through hb annihilation as an example to show 
that it is reasonable to use the two cut-off phase space slicing method in our NLO QCD 
calculations, i.e. the dependence of the NLO QCD predictions on the arbitrary cut-offs 5s and 



5c is indeed very weak, as shown in Ref. 



2q . Here cXother includes the contributions from the 
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Born cross section and the virtual corrections, which are cut-off independent. Both the soft 
plus hard colhnear contributions and the hard non-colhnear contributions depend strongly 
on the cutoffs and, especially for the small cut-offs {6s < 10^^), each is about ten times larger 
than the LO total cross section. However, the two contributions {(Jgoft + <^hard/coii + (^virtual 
and c hard /non- coll) nearly cancel each other completely, especially for the cut-off 6s between 
5 X 10"^ and 10~^, where the final results for aNio are almost independent of the cut-offs 
and very near 7.7 fb. Therefore, we will take 6s = 10~^ and 6c = 5^/50 in the numerical 
calculations below. 

In Figs. IBHOl we give the total cross sections for h^h^, H^H^, A^A^ and h^H^ production, 
respectively, and compare the 66-annihilation contributions with the gfgf-fusion contributions 
m Qj3, wMC. ..e fro. ,.a* a„d .,ua* loop. I„ ca. of p.oduc.o. 
Ref. [iS^ indicated that hb annihilation can be more important than gg fusion for large values 
of tan/9, but it is not so here (see Fig. Efa)), which is due to the fact that we have used a 
much larger decay width for the than the one in Ref. . However, when tan /3 is small 
(< 15), the contributions of hb annihilation still can exceed those of gg fusion (see Fig.E^b)). 
As for H^H^ (Fig- IZ|) and A^A^ (Fig- El) production, hh annihilation is suppressed by a factor 
between 2 and 3 in most of the parameter space compared to gg fusion except for m^o > 400 
GeV and tan/3 = 40, where the contributions of hh annihilation are larger than those of gg 
fusion, but the corresponding cross sections are very small (< 0.1 fb). In the case of H^h? 
production, hh annihilation dominates for m^o > 250 GeV and large values of tan (3. From 
Figs. EHHl we also see that the NLO QCD corrections to the total cross sections for these 
four processes can enhance the LO results significantly for > 0, generally by a few tens of 
percent, while for /i < 0, the corrections are relatively small, and are even negative in some 
parameter regions. 

In Figs. ITUl and ITTl we plot the total cross sections for pp A^h^, A^H^ at the LHC 
as functions of m^o and tan /3, and compare the hh annihilation contributions with gg fu- 
sion and qq annihilation [li|. In the case of A^h^ production, the hh- 
annihilation contributions dominate for large values of tan (3. For example, when tan /? = 40 
and m^o > 250 GeV (see Fig. [TUT a)). the contributions are several times larger than gg 
fusion contributions, and at least two orders of magnitude larger than qq annihilation con- 
tributions. However, for small values of tan/3, gg fusion dominates, as shown in Fig. [TUT b). 
In the case of A^H^ production, after including the NLO QCD corrections, the cross section 
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for bb annihilation is lager than those of the other two mechanisms for large values of tan (3 
and most values of rriAo. Moreover, the 66-annihilation contributions can exceed 100 fb for 
tan/3 = 40 and m^o < 150 GeV, as shown in Fig. ITlT a). From Figs. ITUl and fTH we also 
see that the NLO QCD corrections to the total cross section for these two processes can 
enhance the LO results significantly for /i > 0, generally by a few tens of percent, while for 
/i < 0, the corrections are negative and relatively small. 

Fig. ^1 gives the dependence of the ratio K (defined as the ratio of the NLO total cross 
sections to the LO ones) on m^o for A^h^/A^H^ production through bb annihilation based 
on the results shown in Fig. imT a) and HIT a). We see that in general the ratio K is negative, 
and becomes larger with the increasing m^o. For example, when m^io varies from 120 GeV 
to 500 GeV, the ratio K increases from 0.7 to 0.95 for A^h^ production (Fig. Efa)), and 
from 0.75 to 0.88 for A^H^ production fFig. [T^ b)). The contributions to the ratio K come 
from three IR finite parts: the LO total cross sections, the pure QCD corrections and the 
SUSY-QCD virtual corrections, the latter two of which are also shown in the figure. 

The main parts of bb annihilation contributions for large values of tan (3 originate from the 
Yukawa coupling b — b — (f), so the results are sensitive to the bottom quark mass. In Fig. 
we show the effects of the choices of on the total cross section for A^H^ production 
through bb annihilation, assuming tan/3 = 40, mi/2 = 170 GeV, Aq = 200 GeV and /i < 0. 
We considered three different bottom quark masses, (i) MS bottom quark mass at the scale 
of the MS mass, i.e. mbirrib) = 4.25 GeV, (ii) QCD running bottom quark mass rriblfir), 
and (iii) QCD running plus SUSY improved bottom quark running mass. They can have 
very different values. For example, in this figure, when m^o = 250 GeV, they are 4.25 GeV, 
2.69 GeV and 3.17 GeV, which leads to the LO total cross sections being 12.5 fb, 5 fb and 
9.8 fb, respectively. Moreover, from Fig. ^[ we see that the NLO QCD corrections can be 
very large in the cases of mb{mb) and QCD running mf,. For this reason we use the SUSY 
improved bottom quark running mass, which improves the convergence of the perturbation 
calculations, especially for large values of tan /5, as shown in Fig. El 

Fig. El gives the dependence of the total cross section for A^H^ production through hb 
annihilation at the LHC on the renormalization/factorization scale for /x^ = A*/- In the case 
of /i > 0, the scale dependence of both the LO and the NLO total cross sections is relatively 
weak. And for < 0, the scale dependence of the total cross sections is reduced when going 
from LO to NLO. For example, the cross sections vary by ±20% at LO but by ±13% at 
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NLO in the region 0.5 < fif/mg^v < 2.0. 

Since another source of uncertainty arises from the different choice of PDFs, in Fig. 
we show the total cross sections for production through bb annihilation as functions 

of rriAO for three different PDFs. We first use the 41 CTEQ6.1 PDF sets 0] to estimate 
the uncertainty in the LO total cross sections. The LO results using the CTEQ6M PDFs lie 
between the maximum and the minimum. The NLO total cross sections are then calculated 
using three different PDF sets, one of which is CTEQ6M, and the other two are the ones 
that gave the maximum and minimum LO uncertainties. Observe that in this case the 
uncertainty arising from the choice of PDFs increases with the increasing m^o. Moreover, 
the dependence of the total cross sections on PDFs is not decreased from LO to NLO. 

In Fig. E[ we display differential cross sections as the functions of the transverse momen- 
tum pt of the and the invariant mass Mj^of^o, which are given by Eqs. ()3.48p and ()3.52p . 
respectively, for the A^H^ production through bb annihilation. In the case of /i < 0, we find 
that the NLO QCD corrections reduce the LO differential cross sections except for low px, 
while in the case of /i > 0, the corrections always enhance the LO results. 

In conclusion, we have calculated the complete NLO inclusive total cross sections for pair 
production of neutral Higgs bosons through bb annihilation in the MSSM at the LHC In 
our calculations, we used both the DREG scheme and the DRED scheme and found that 
the NLO total cross sections in the above two schemes are the same. Our results show 
that the bb annihilation contributions can exceed those of gg fusion and qq annihilation for 
hPH^, A^h*^ and A^H^ production when tan/5 is large. For /i > the NLO corrections 
enhance the LO total cross sections significantly, and can reach a few tens percent, while 
for /i < the corrections are relatively small and negative in most of parameter space. 
Moreover, the NLO QCD corrections reduce the dependence of these total cross sections on 
the renormalization/factorization scale, especially for /i < 0. We also used the CTEQ6.1 
PDF sets to estimate the uncertainty in both the LO and NLO total cross sections, and 
found that the uncertainty arising from the choice of PDFs increases with the increasing 
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APPENDIX A 



In this appendix, we give the relevant Feynman couplings j41|. 

1. Hi^b-b: igYbia^PL + a*PR) 



ai 



; COS a, 



02 



sm a, 



as 



1 1 

-j=smP, 04 = -^ cos/?. 



v/2 ' ^ V2 
where a is the mixing angle in the CP-even neutral Higgs boson sector. 
2. Hj-Z-A^: gZfipAo - pH,Y/i2cw) 



0. 



= — sin(/3 — a), Z2 =cos(/? — a), ^3 

Here we define the outgoing four-momenta of and Hj positive. 

3. Hk- Hi - Hj : igmzCkijl {2cw) 

Cm = — 3 cos 2a; cos(a; + P), C112 = 2sin2acos(a + P) + sin(a + /3) cos2q;, 
C122 = —2 sin 2a sin(a + P) + cos{a + P) cos 2a, C133 = cos 2/3 cos(a + /?), 
C134 = sin 2/? cos(a + /?), C144 = — cos 2/5 cos(a + /5), 
C222 = —3 cos 2a sin(a + /?), C233 = — cos 2/3 sin(a + /?), 
C244 = cos 2/5 sin(a + /?), C234 = — sin 2/? sin(a + /?). 

The indexes i,j and k of Ckij are symmetric, and other coefficients are zero. 

4. Hk-bi- bm-. igGt = tg[R'G'^{R'f]im 

-— cos(a + P)CbL — V^rrihYbCosa — j=Yb(Ai,cosa — /isina) 

-^Yb{Ati cos a — /isina) ^ cos(a -|- P)CbR — v^mfeYfe cos a 



Cw 



I 



G' 



— sin(a + P)CbL + "\/2"^f,Yf, sin a -j=Yb{Ab sin a + /i cos a) 

-^Yi,{AbSma + //cos a) sin(a + P)GbR + v^mfeYj, sin 



a 



20 



rUb 



—Af, tan [3 — 11 

'^^w [ Abtan(3 + fi 



G' 



rrib 



2m 



w 



( 



Ab — ^ tan l3 

—Ab + fitan/3 



5. Z^-k- L: -tgZUp-b^ +Pir/cw 



where and pj,^ are the four-momenta of hm and hm in direction of the charge flow. 
6. Hi - H, -hi- K: tg^Gl^ = tg^[R^G^^ {R'Y]im 



All 



A22 



A33 



A44 



Al2 



A34 



1 / COS 2a m^cos^a \ r t- 

2l^|-^^^ + mf,cos2/?j^'^^-^^2 



1 / cos 2a 



1 / COS 2a 

2 ^ 



ChT, — 



C 



w 



1 /cos 2/5 



r2 



6L 



I'W 

ml sin^ a 
cos2/5^ 

m? tan^ /5\ 



cos a 
m|(/ cos^ 

,2 



I20m2i 



1 / COS 2a sm' a > 

-(--fe/? H ^ TT; I Ol20m2, 



c 



w 



1 /cos 2/5 



mt 



1 /cos2/5 

CnCml — - I 5 ^bR + 



cos^ /5 
ml tan^ /3\ 



r2 



m 



5«2'5m2, 



^/l^ml + 



1 /cos 2/? 



sin 2a 
2 

sin 2/? 



+ 



r2 



CbR — 



2m^cos2/5 



bL 



+ 



2m^ cos^ /5^ 



'G 



bR 



mb 
mt 



SioS. 



720m2 5 



r2 



bR 



2m^ cos^ P ^ 
^1 

2m^r cos^ /5 



720m2 



i20m2 



The indexes i and j of G]l^ are symmetric, and other coefficients are zero. 
APPENDIX B 

In this appendix, we collect the explicit expressions for the nonzero form factors in 
Eq. ()3.6p . For simplicity, we introduce the following abbreviations for the Passarino-Veltman 
integrals, which are defined as in Ref. j2J] except that we take internal masses squared as 
arguments: 



B^'' =Bo{a,0,0), 

,2 2 



5?"^ = Bo{a,ml,m- J, 
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Cp%) = Cp(g)(0, 0,5,0,0,0), 

Cpig) = Cp(<7)(0, a, 0, ml ml^mlj, 

C^pt;^ = Cp(.)(0,<.,«, 0,0,0), 



^d{ar) rf 1 1\ 2 2 2 2 \ 



= D 



D 



b(a) 



Pig) 



a, m 



2 ,s, 0,0, m^^, 0,0, 0,0), 



0(a 7-1 /n 2 n2 222 2\ 



Many functions above contain soft and collinear singularities, but all the Passarino-Veltman 
integrals can be reduced to the scalar functions Bq, Cq and Dq. Here we present the explicit 
expressions for Cq and Dq, which contain the singularities, and were used in our calculations: 



27r2 



a — m 



Hr 



lln 



-a 



m 



m 



Hr 



2 \-a 



— 

sa 
-Li2 



12 fniH.mH^ 
^ + - In ' 



e2 e 



-a 



+ 



27r^ 



— sa 



s-m'\j^){m]j^ - a) 
—sa 



3 

+ Li2 



-a 



TT 
.2 



— I Li 



Li2 



s — m 



m 



Li2 



2 L(s-m|^J(m|^^ -a 
1 —2 
2 



+ -ln^ 



—a 



+ In 



In 



s + a — m 



—a 



a^ 



1 

4 



1 
2 



+- In ( — ^ I In ' 



m 



1 



a" 
m 



2 Kttijj, — a 



where we define = (47r/i^/s)T(l — e)/r(l — 2e), and note the explicit expression for Dq 



a{a) 



is in agreement with the one in Ref. 



421. 



There are the following relations between the form factors: 



/2,4 ~ fl,3i^i,. 



Thus we will only present the explicit expressions of and f^. Corresponding to diagrams 
(a)-(m) in Fig. |2l the form factors are 

n = E ^^^{na.[s(Co" + 2C1 + (2 - e)C1,) + 4(6 - 2)CSq + 1] 
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2 

l,m=l 

l,m=l ^W\^ "^ZJ 

fl = .2^. '^2 ^ ^^^[^(2^0" + 4Ci« + 2(1 - 6)C,»,) - 4(1 - 6)C«o - 1] 

2 

i,m=l 

2 

/i + /r = - H '^Yb{ajG]jnRiiR^^iC2^^'^ + aiG'^;i?|'2-R^2C'f 

Lm=l 



r=t,] 

+Clt^ + C^'^*''^) + 4(2 - e)Cof - m^^(Co^(*^) + cf'^ + (2 - e)Cit'') + Cf 

2 

~1] + X! "^9("i^/m-^«2-RmlC'o ' + aj<^m«-^ml-^i2C'o )}, 

i,m=l 

+ // - - E 2n(a,G'^,4^L2^2^"'^ + a.GL4^LiG2^"^'^), 

Z,m=l 

+C2t''^ + Cl^""'^) + 4(2 - e)(7ot^) - m^,^ (C^^^'"') + C^^^^^ + (2 - e)Clt''^ + C^'^"'')) 

2 

~1] + E ''^oVh^lmRnRml^O + %G^mi-^a-^m2C'o ]}, 
(,m=l 

/f + Z!^ = - E 2n^a,a,44^(t5S^") + X^SS^*)), 

/I + /a'^ = n^{^(to,a*So"H - ^.a.a*<)) + E(4)^[^(K " 

fi = 2y,2a,a*[Co^(**) + Cq^^*^') - (1 + e)Cl^'^ - sD^^ + {u - te)Df\ 

i,m,n=l 
2 

l,m,n=l 

= -2Y^aja*[Cl^'''^ + Cq"^"^') - (1 + e)C'^^'^ - sDl^""^ + {t - ue)Df% 

2 

<3 r>b r>b r>b{u) 



fi — ~ E 2mgG'^„G^;it!^24i-^o > 

i,m,n=l 
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Js — ^'-^mn'-'n/-'T-/l-'T-ml-'^3 5 



/,r7i,n=l 
2 

/l — '^"^5'-^Zm-'T^/2-'T^ml'-^0 • 

Lm=l 



APPENDIX C 

In this Appendix, we collect the explicit expressions for the amplitudes squared for the 
radiation of a real gluon.The results for massless b{b) emission can be obtained by crossing 
symmetry. Since these expressions are only used for the hard non-collinear parts of the real 
corrections in Eq. (j3.18| ) and Eq. (|3.33|) . they have no singularities and can be calculated in 
n = 4 dimensions. 

For simplicity, we define the following invariants: 

S3 = iP2 + Pzf, Sa = {pi + psf, S5 = {pi + P2)^ 

t={Ph-p2f, t'={pb-p^f, U={pa-P2f, 

U=iPa-P3f, %=(Pb-Pl)^ Uj={pa-Plf, 

. Then for real gluon emission we find 



with 



^HH)r 2m|n2 ^ is + t'-u%s~t' + u') 



\Mij I — 2 X! + '^*k'^l)CkijClij Tn T 

k,l=l '^'^^Hk^H, 

'\ = ' 2 \2+, / l-^4M6(g + t} + S4U7{S + u)+ U(iUj{2s + t + M ) - t ».y 

— Mg-u' — mjj_ [t'(s4 + t' + 2uq — uj) + u\u + S4 — Me + 2M7) + s(2s + 2s4 + 2t' 

+3% + 3u7 + 2u') + m^^(4s + t' + m')]}, 
-2 ^5,3n'^f . „ ... 2 



Wif^^ = n ^ % . {CbRaja* + CbLaia*){-£—{s + t' + uq - )[mi (2s + 2s4 + t' 
+2(m6 + 2u7 + u) - 4:771^,) - (2-S4 + t' + 2uq)u7] 

— - — [ss4 + ^7(454 + t') — uqu — m}j^{s + 4s4 + t' + Auj + 3m' — 4m^.)] 

'U7 

^ ;(m7 - m]j^Wsi + 2m7 + 3m' - Am]j;){m]j^ - ue) - {s + t'){u' - 2mjjJ] 



tU'jU 
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H YIuq — 2uj — u)(uqU — t'uY + (t' — u )'mjT ) + s(3u + m\ (4s + At' 

+Q{uq + M7) — AuqUj — 8m^.) + Sii^u' — 2uq — 2^7 + 4m^.))] 
-^[53(5 + 4t) - t'u + tu' - m]j^{s + 4s3 + 4t + 3t' + u' - 4m^^.)] 

+ ^[s(s3(t' - 2t - 2m) - Atu) + (2t + t' - 2u){tu - t m) + ((2t + t' 
-2u){t' - u) + s(4s + 4s3 + 6t + 3t' + 6m + Au - 8m^,))]}, 



|Mj^j- '^l = leyj, |ai| |aj| {—7^(5 + t' + uq - mfj.)[u-j{si + Me) - m^j.^s + S4 + Me + 2^7 
+u' - 2m]j)] - ^ [(Me - - M')(MeM' - t'M7 + m]j it' - u')) 

tt'UjU' 

— s(2m6M7 + S4(m6 + M7 — m')) + sm^. (2s + 2s4 + 2t' + 3uq + 3m7 + m' — 4m^J] 
2 

+ ^^^^(^7 - "^hJ[(s + t' + 2m6 - 2m^J(m^^ - m') + (m^^ - Me) (54 + M7)] 
+:^[m7((s + t')(s4 + ur) - Uqu') - m]j.{t\si + u' - mjjj + ujis + 2t' - u'))] 

6 

+^[^3^ - mjj^iss + t + t' - m^^)] + J^isiUj - m^^(s4 + M7 + m' - m^J]}, 

\Mlp\ = lQY^alaf{-^\t'm - u^u - (34 - m\^{s + 2m6) + m^^(2s4 + t' + m' - 2m^J] 

2 2 

+ 7 + ^ (tMe - m|^^mi.)[M(M7 - m^J + M7(m - m^J] 

-[(s4 — m'jjj{s + 2uj) + t'M7 — u^u — m\.{2si + t' + u — 2m^.)] 



t'MM7 

2 

~ t'uuGU7 ^^^ - ^nM^m ~ "7)(s4 + + Me) + + 2m7 + u' - 2mj^J] 

+ Tr~ ;[(^6 - Uj){uqU - t'u-j) - s(2m6M7 + S4(Me + M7)) + m%.{{uQ - M7)(t' - u ) 

2 

+s{2s + 2s4 + 2t' + 3m6 + 3m7 + 2m' - 4m^J)] + — ; (t - mjj,){uuj - m]jjn]j) 

+ — ^[(t - M)(tM' - t'u) - s{2tu + ssit + m)) + mjr {{t - u){t' - u') + s{2s 
tt'uu ^ 

+2s3 + 3t + 2t' + 3m + 2m' - 4m|,))]}. 
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FIG. 1: Leading order Feynman diagrams for bb — > HiHj. 
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FIG. 3: Feynman diagrams for the real gluon emission contributions. 




FIG. 4: Feynman diagrams for the massless bottom quark emission contributions. 
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FIG. 5: Dependence of the total cross sections for A^H^ production through bb annihilation at the 
LHC on the cut-off Sg, assuming m^o = 250 GeV, tan/3 = 40, mi/2 = 170 GeV, Aq = 200 GeV, 
fi <0 and Sc = 6s/50. 
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FIG. 6: Total cross sections for the h^h^ production at the LHC (a) as functions of m^o for 
tan/3 = 40, /X < (starting from m^o = 120 GeV) and tan /3 = 40, /x > (starting from m^o = 210 
GeV), and (b) as functions of tan/3 for m^o = 250 GeV, assuming mi/2 = 170 GeV and Aq = 200 
GeV. 
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FIG. 7: Total cross sections for the H^H^ production at the LHC (a) as functions of m^o for 
tan /3 = AO, fi < (starting from m^o = 120 GeV) and tan /? = 40, > (starting from m^o = 210 
GcV), and (b) as functions of tan/? for m^o = 250 GeV, assuming mi/2 = 170 GeV and = 200 
GeV. 
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FIG. 8: Total cross sections for the A^A'^ production at the LHC (a) as functions of m^o for 
tan/3 = 40, < (starting from m^o = 120 GeV) and tan /3 = 40, /x > (starting from m^o = 210 
GeV), and (b) as functions of tan/? for m^o = 250 GeV, assuming mi/2 = 170 GeV and = 200 
GeV. 
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FIG. 9: Total cross sections for the production at the LHC (a) as functions of m^o for 

tan/3 = 40, /X < (starting from m^o = 120 GeV) and tan /3 = 40, /x > (starting from m^o = 210 
GeV), and (b) as functions of tan/3 for m^o = 250 GeV, assuming mi/2 = 170 GeV and Ao = 200 
GeV. 
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FIG. 11: Total cross sections for the A^H^ production at the LHC (a) as functions of m^o for 
tan/? = 40, /X < (starting from m^o = 120 GeV) and tan /3 = 40, /x > (starting from m^o = 210 
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FIG. 13: Dependence of the total cross sections for the A^H^ production through bb annihilation 
at the LHC on m^o and mf,, assuming tan/3 = 40, mi/2 = 170 GeV, Aq = 200 GeV and fi < 0. 
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FIG. 14: Dependence of the total cross sections for the A'^H^ production through bb annihilation at 
the LHC on renormalization/factorization scale, assuming tan/? = 40, mi/2 = 170 GeV, = 200 
GeV and /x < 0. 
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